Introducción

En este documento se aplica la regresión Poisson para predecir el número de vehículos registrados diariamente en Colombia en el Registro Único de Tránsito (RUNT).

Los datos son obtenidos del RUNT a través de la compra de información y se muestran aquí solo con fines educativos.

Hipótesis de trabajo

El número de unidades registradas en el RUNT está relacionada con la actividad comercial y de las oficinas de tránsito, que operan mayoritariamente entre lunes y viernes. Además, las fechas especiales pueden tener un impacto.

Lectura de los datos

La base de datos original tiene solo dos campos: fecha y unidades registradas en el RUNT en esa fecha. La base de datos se aumenta utilizando variables indicadoras para los días de la semana, los meses del año y las fechas especiales. Esta es la base que se lee a continuación.

Lectura de los datos

Antes de leer los datos de un archivo externo, y si estos no ocupan un volumen importante, es conveniente visualizarlos para saber qué tipos de datos se tienen y cómo leerlos mejor.

En este caso se leerán los datos con el paquete readr y la función read_delim():

library(readr)
datos_ventas_autos <- read_delim("datos_ventas_autos.csv", 
     ";", escape_double = FALSE, trim_ws = TRUE)
head(datos_ventas_autos) # muestra una vista de los primeros 10 registros

La función header() muestra una vista de los primeros 6 registros.

Análisis descriptivo

El análisis descriptivo busca visualizar comportamientos de los datos como dispersión, simetría o asimetría, relaciones entre variables, etc., mediante el cálculo de unas variables de resumen, llamadas estadísticos, y la representación gráfica de los datos.

Para guiar el análisis de la información resulta útil representar gráficamente los datos o una muestra de los mismos. En este caso, como la variable a predecir es el número de unidades registradas diariamente en el RUNT y esta variable depende del tiempo, se desea representarla como una serie temporal.

Sin embargo, la variable Fecha es una variable tipo caracter:

class(datos_ventas_autos$Fecha)
## [1] "character"

Para representar el número de unidades registradas en el RUNT en el tiempo es útil cambiar el tipo de variable a caracter. El siguiente código hace esto usando la función as.Date():

datos_ventas_autos$Fecha<-as.Date(datos_ventas_autos$Fecha,"%d/%m/%Y")

Ahora se procede a graficar los datos para darse una idea de algunos comportamientos:

library(plotly)
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~Fecha,
         y = ~und,
         type = "scatter" ,mode = "lines",
         line=list(width=1,color='rgb(205, 12, 24)'))%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Día"),
         yaxis=list(title="Unidades"))

Ahora veamos el comportamiento anual. Creemos la columna “year” (año):

datos_ventas_autos$year<-format(datos_ventas_autos$Fecha,"%Y")
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~Fecha,
         y = ~und,
         type = "scatter" ,mode = "lines",
         split = ~year,
         line=list(width=1))%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Día"),
         yaxis=list(title="Unidades"))

Otra forma de analizar los resultados por año es con la función aggregate(). Utilicemos esta función para obtener el promedio diario para cada año de registros de vehículos:

aggregate(und~year,data=datos_ventas_autos,FUN=mean)

Ahora obtengamos el promedio diario para cada mes y cada año:

datos_ventas_autos$mes<-strftime(datos_ventas_autos$Fecha, format = "%B")
datos_ventas_autos$mes<-factor(datos_ventas_autos$mes,levels = month.name)
aggregate(und~year*mes,data=datos_ventas_autos,FUN=mean)

Veamos este promedio representado en un gráfico:

aggregate(und~year*mes,data=datos_ventas_autos,FUN=mean)%>%
  plot_ly(x = ~mes,
         y = ~und,
         type = "scatter" ,mode = "lines",
         split = ~year,
         line=list(width=1))%>%
  layout(title='Promedio diario mensual de registros el en RUNT',
         xaxis=list(title="Día"),
         yaxis=list(title="Unidades"))

Ahora obtengamos el día del año para comparar todos los años:

datos_ventas_autos$diaano<-strftime(datos_ventas_autos$Fecha, format = "%j")

Obtengamos la gráfica del registro de vehículos diarios para cada año:

library(plotly)
plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~diaano,
         y = ~und,
         type = "scatter" ,mode = "lines",
         split = ~year,
         line=list(width=1))%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Día"),
         yaxis=list(title="Unidades"))

Ahora veamos el comportamiento mensual. Primero se extrae la información del mes y se expresa como una variable categórica ordenada:

datos_ventas_autos$diames<-format(datos_ventas_autos$Fecha,"%d")

La siguiente función “automatiza” los pasos para generar un gráfico para cada mes:

grafico<-function(month,data,sl=FALSE){
  plot_ly (data=subset(data,subset = (Fecha<="2017-12-31" & mes==month)),
         x = ~diames,
         y = ~und,
         type = "scatter" ,mode = "lines",
         split = ~year,
         line=list(width=1),
         legendgroup=~year,
         showlegend = sl)%>%
    layout(title='Registros el en RUNT',
         xaxis=list(title="Día"),
         yaxis=list(title="Unidades"))->p
  return(p)
}

Apliquemos esta función a todos los meses del año:

meses<-unique(datos_ventas_autos$mes)
gr1<-grafico(meses[1],data=datos_ventas_autos,sl=FALSE)
gr2<-grafico(meses[2],data=datos_ventas_autos,sl=FALSE)
gr3<-grafico(meses[3],data=datos_ventas_autos,sl=FALSE)
gr4<-grafico(meses[4],data=datos_ventas_autos,sl=FALSE)
gr5<-grafico(meses[5],data=datos_ventas_autos,sl=FALSE)
gr6<-grafico(meses[6],data=datos_ventas_autos,sl=FALSE)
gr7<-grafico(meses[7],data=datos_ventas_autos,sl=FALSE)
gr8<-grafico(meses[8],data=datos_ventas_autos,sl=FALSE)
gr9<-grafico(meses[9],data=datos_ventas_autos,sl=FALSE)
gr10<-grafico(meses[10],data=datos_ventas_autos,sl=FALSE)
gr11<-grafico(meses[11],data=datos_ventas_autos,sl=FALSE)
gr12<-grafico(meses[12],data=datos_ventas_autos,sl=TRUE)
subplot(gr1,gr2,gr3,gr4,gr5,gr6,gr7,gr8,gr9,gr10,gr11,gr12,nrows = 4,shareX = TRUE,
        shareY = TRUE,which_layout=1)

Una alternativa, más simple pero con menos interacción, para la gráfica anterior se obtiene usando la función xyplot():

library(lattice)
xyplot(und~diames|mes,groups = year,data=datos_ventas_autos,
       subset = (Fecha<="2017-12-31"),type=c("l","g"),
       auto.key = list(lines = TRUE, space = "bottom",columns=7,title="Año",cex=0.5),
       layout=c(3,4))

Ahora utilicemos el diagrama de caja y bigotes para explorar relaciones:

plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~year,
         y = ~und,
         type = "box")%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Año"),
         yaxis=list(title="Unidades"))

Ahora, hacemos algo similar a lo anterior para los meses del año y los días de la semana:

datos_ventas_autos$dia_semana<-weekdays(datos_ventas_autos$Fecha)
datos_ventas_autos$dia_semana<-ordered(datos_ventas_autos$dia_semana,levels=c( "Tuesday", "Wednesday", "Thursday", 
"Friday", "Saturday", "Sunday","Monday"))

Veamos el diagrama de caja y bigotes para cada mes:

plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~mes,
         y = ~und,
         type = "box")%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Mes"),
         yaxis=list(title="Unidades"))

Veamos el diagrama de caja y bigotes para cada día de la semana:

plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~dia_semana,
         y = ~und,
         type = "box")%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Día de la semana"),
         yaxis=list(title="Unidades"))

¿Cuál será el efecto de los días feriados?

plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~Feriado,
         y = ~und,
         type = "box")%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Día feriado (1 es Si)"),
         yaxis=list(title="Unidades"))

¿Cuál será el efecto del día del amor y la amistad?

plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~AmoryAmistad,
         y = ~und,
         type = "box")%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Amor y amistad (1 es Si)"),
         yaxis=list(title="Unidades"))

¿Cuál será el efecto de la Semana Santa?

plot_ly (data=subset(datos_ventas_autos,subset = (Fecha<="2017-12-31")),
         x = ~`Semana Santa`,
         y = ~und,
         type = "box")%>%
  layout(title='Registros el en RUNT',
         xaxis=list(title="Semana santa (1 es Si)"),
         yaxis=list(title="Unidades"))

Modelamiento

Se ajustarán modelos con la información disponible hasta el 31 de diciembre de 2016 y se utilizará el año 2017 para validar el modelo:

Regresión lineal:

Utilizamos la función lm():

modelo_lm<-lm(und~dia_semana+mes+Feriado,data=datos_ventas_autos,subset = (Fecha<="2016-12-31"))

Veamos el resultado:

summary(modelo_lm)
## 
## Call:
## lm(formula = und ~ dia_semana + mes + Feriado, data = datos_ventas_autos, 
##     subset = (Fecha <= "2016-12-31"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1077.39  -116.40   -20.32   108.56  2096.80 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   682.415     21.479  31.771  < 2e-16 ***
## dia_semana.L -698.583     16.790 -41.608  < 2e-16 ***
## dia_semana.Q  -28.194     16.674  -1.691 0.091025 .  
## dia_semana.C  729.146     16.603  43.915  < 2e-16 ***
## dia_semana^4  491.497     16.509  29.772  < 2e-16 ***
## dia_semana^5    8.992     16.495   0.545 0.585701    
## dia_semana^6 -225.718     16.492 -13.687  < 2e-16 ***
## mesFebruary   117.821     31.003   3.800 0.000149 ***
## mesMarch      174.696     30.267   5.772 9.21e-09 ***
## mesApril      168.802     30.511   5.533 3.62e-08 ***
## mesMay        147.077     30.258   4.861 1.27e-06 ***
## mesJune       145.742     30.511   4.777 1.93e-06 ***
## mesJuly       153.337     30.261   5.067 4.45e-07 ***
## mesAugust     167.486     30.258   5.535 3.56e-08 ***
## mesSeptember  165.400     30.570   5.411 7.12e-08 ***
## mesOctober    150.525     30.273   4.972 7.25e-07 ***
## mesNovember   171.269     30.509   5.614 2.29e-08 ***
## mesDecember   464.365     30.259  15.346  < 2e-16 ***
## Feriado      -819.039     28.770 -28.469  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.4 on 1808 degrees of freedom
## Multiple R-squared:  0.7696, Adjusted R-squared:  0.7673 
## F-statistic: 335.5 on 18 and 1808 DF,  p-value: < 2.2e-16

Cálculo del \(R^2\) manualmente:

y0_tr<-mean(datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"])
r0_tr<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y0_tr
R0_tr<-mean(r0_tr^2)
y_pred_tr_lm<-predict(modelo_lm)
r_tr_lm<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y_pred_tr_lm
R_tr_lm<-mean(r_tr_lm^2)
R2_lm<-1-R_tr_lm/R0_tr
print(R2_lm)
## [1] 0.7695787

Ajuste de un modelo Poisson:

modelo_glm<-glm(und~dia_semana+mes+Feriado,data=datos_ventas_autos,subset = (Fecha<="2016-12-31"),family = "poisson")

Veamos el resultado:

summary(modelo_glm)
## 
## Call:
## glm(formula = und ~ dia_semana + mes + Feriado, family = "poisson", 
##     data = datos_ventas_autos, subset = (Fecha <= "2016-12-31"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -42.581   -3.415   -1.120    2.933   45.878  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   5.795637   0.005850  990.71   <2e-16 ***
## dia_semana.L -2.559408   0.013146 -194.69   <2e-16 ***
## dia_semana.Q  0.101049   0.002206   45.81   <2e-16 ***
## dia_semana.C  2.800453   0.014192  197.32   <2e-16 ***
## dia_semana^4  3.120679   0.019455  160.40   <2e-16 ***
## dia_semana^5  1.858125   0.015167  122.51   <2e-16 ***
## dia_semana^6  0.524532   0.007138   73.49   <2e-16 ***
## mesFebruary   0.168204   0.004353   38.64   <2e-16 ***
## mesMarch      0.219482   0.004304   51.00   <2e-16 ***
## mesApril      0.217067   0.004310   50.36   <2e-16 ***
## mesMay        0.208480   0.004292   48.57   <2e-16 ***
## mesJune       0.203179   0.004353   46.67   <2e-16 ***
## mesJuly       0.207886   0.004270   48.68   <2e-16 ***
## mesAugust     0.234499   0.004276   54.84   <2e-16 ***
## mesSeptember  0.226172   0.004252   53.19   <2e-16 ***
## mesOctober    0.208291   0.004251   48.99   <2e-16 ***
## mesNovember   0.235782   0.004310   54.70   <2e-16 ***
## mesDecember   0.545568   0.004011  136.03   <2e-16 ***
## Feriado      -5.662272   0.059671  -94.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 913454  on 1826  degrees of freedom
## Residual deviance: 108666  on 1808  degrees of freedom
## AIC: 121839
## 
## Number of Fisher Scoring iterations: 8

Cálculo del Pseudo \(R^2\) para el modelo Poisson:

Obtengamos el (pseudo) \(R^2\) para el modelo Poisson:

y_pred_tr_glm<-predict(modelo_glm,type="response")
r_tr_glm<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y_pred_tr_glm
R_tr_glm<-mean(r_tr_glm^2)
R2_tr_glm<-1-R_tr_glm/R0_tr

Veamos el resultado:

print(R2_tr_glm)
## [1] 0.807181

Análisis de predichos versus observados

Creación de un dataframe para usar con plotly:

resultados_lm_glm<-data.frame(Fecha=  datos_ventas_autos$Fecha[datos_ventas_autos$Fecha<="2016-12-31"],                                    und=datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"],
                              pred_lm=y_pred_tr_lm,
                              pred_glm=y_pred_tr_glm,
                            res_lm=r_tr_lm,
                            res_glm=r_tr_glm)

Predichos y observados:

plot_ly (data=resultados_lm_glm,
         x = ~Fecha,
         y = ~und,
         type = "scatter" ,mode = "lines",
         name='Real',
         line=list(width=1,color='rgb(205, 12, 24)'))%>%
  add_trace(y= ~pred_lm,
            name='Modelo lineal general',
            line=list(width=1,color='rgb(22, 96, 167)'))%>%
  add_trace(y= ~pred_glm,
            name='Modelo Poisson',
            line=list(width=1,color='rgb(255, 51, 0)'))%>%
  layout(title='Registros RUNT',
         xaxis=list(title="Fecha"),
         yaxis=list(title="Unidades"),
         legend = list(x = 0.75, y = 0.9))

Predichos vs observados:

plot_ly (data=resultados_lm_glm,
         x = ~und,
         y = ~pred_lm,
         text = ~Fecha,
         type = "scatter" ,mode="markers",
         name='Modelo lineal general',
         marker=list(size=3,color='rgb(22, 96, 167)'))%>%
  add_trace(y= ~pred_glm,
            text = ~Fecha,
            name='Modelo lineal Poisson',
            marker=list(size=3,color='rgb(255, 51, 0)'))%>%
  add_trace(x=c(-500:3500),y=c(-500:3500),
            mode="lines",text=rep(NA,4001),
            name="Identidad")%>%
  layout(title='Registros RUNT',
         xaxis=list(title="Observados"),
         yaxis=list(title="Predichos"),
         legend = list(x = 0.75, y = 0.9))
## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...

Residuales vs observados:

plot_ly (data=resultados_lm_glm,
         x = ~und,
         y = ~res_lm,
         type = "scatter" ,mode="markers",
         name='Modelo lineal general',
         text = ~Fecha,
         marker=list(size=3,color='rgb(22, 96, 167)'))%>%
  add_trace(y= ~res_glm,
            name='Modelo Poisson',
            text = ~Fecha,
            marker=list(size=3,color='rgb(255, 51, 0)'))%>%
  layout(title='Registros RUNT',
         xaxis=list(title="Observados"),
         yaxis=list(title="Residuales"),
         legend = list(x = 0.75, y = 0.9))

Predicción del 2017

Cálculo de las predicciones:

datos_val<-subset(datos_ventas_autos,subset=(Fecha>="2017-01-01" & Fecha<="2017-12-31"))
y_pred_vl_lm<-predict(modelo_lm,newdata = datos_val)
y_pred_vl_glm<-predict(modelo_glm,type="response",newdata = datos_val)
r_vl_lm<-datos_val$und-y_pred_vl_lm
r_vl_glm<-datos_val$und-y_pred_vl_glm
r0_vl<-datos_val$und-y0_tr

Cálculo de los \(R^2\) de predicción:

R_vl_lm<-mean(r_vl_lm^2)
R_vl_glm<-mean(r_vl_glm^2)
R0_vl<-mean(r0_vl^2)
R2_vl_lm<-1-R_vl_lm/R0_vl
R2_vl_glm<-1-R_vl_glm/R0_vl
print(R2_vl_lm)
## [1] 0.5184184
print(R2_vl_glm)
## [1] 0.5511457

Gráfico de los residuales en validación:

resultados_lm_glm_val<-data.frame(Fecha=datos_val$Fecha,                            und=datos_val$und,
                              pred_lm=y_pred_vl_lm,
                              pred_glm=y_pred_vl_glm,
                            res_lm=r_vl_lm,
                            res_glm=r_vl_glm)
plot_ly (data=resultados_lm_glm_val,
         x = ~Fecha,
         y = ~und,
         type = "scatter" ,mode = "lines",
         name='Real',
         line=list(width=1,color='rgb(205, 12, 24)'))%>%
  add_trace(y= ~pred_lm,
            name='Modelo lineal general',
            line=list(width=1,color='rgb(22, 96, 167)'))%>%
  add_trace(y= ~pred_glm,
            name='Modelo Poisson',
            line=list(width=1,color='rgb(255, 51, 0)'))%>%
  layout(title='Registros RUNT (Validación)',
         xaxis=list(title="Fecha"),
         yaxis=list(title="Unidades"),
         legend = list(x = 0.75, y = 0.9))

¿Qué es lo que el modelo no ha visto?

Días hábiles

Creación de un vector que indica si la fecha es el último día hábil del mes. Los días hábiles son aquellos que no sean ni sábado ni domingo ni festivo:

dias_habiles<-ifelse(!(datos_ventas_autos$dia_semana %in% c("Saturday","Monday")) | 
                       datos_ventas_autos$Feriado==0,1,0)
con_mod<-datos_ventas_autos$Consecutivo*dias_habiles
ultimo_dia_habil<-aggregate(con_mod~datos_ventas_autos$year*datos_ventas_autos$mes,FUN=max)
datos_ventas_autos$ultimo_dia_habil<-(datos_ventas_autos$Consecutivo %in% ultimo_dia_habil$con_mod)*1

Ahora incluyamos esta variable en el modelo:

modelo_glm<-glm(und~dia_semana+mes+Feriado+ultimo_dia_habil,
                data=datos_ventas_autos,
                subset = (Fecha<="2016-12-31"),
                family = "poisson")
summary(modelo_glm)
## 
## Call:
## glm(formula = und ~ dia_semana + mes + Feriado + ultimo_dia_habil, 
##     family = "poisson", data = datos_ventas_autos, subset = (Fecha <= 
##         "2016-12-31"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -51.846   -3.189   -0.814    3.231   45.198  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       5.782265   0.005853  987.99   <2e-16 ***
## dia_semana.L     -2.561141   0.013146 -194.83   <2e-16 ***
## dia_semana.Q      0.100992   0.002206   45.78   <2e-16 ***
## dia_semana.C      2.799723   0.014192  197.27   <2e-16 ***
## dia_semana^4      3.119964   0.019455  160.36   <2e-16 ***
## dia_semana^5      1.857067   0.015167  122.44   <2e-16 ***
## dia_semana^6      0.521896   0.007138   73.12   <2e-16 ***
## mesFebruary       0.164640   0.004354   37.81   <2e-16 ***
## mesMarch          0.220647   0.004304   51.27   <2e-16 ***
## mesApril          0.214181   0.004311   49.69   <2e-16 ***
## mesMay            0.208480   0.004292   48.57   <2e-16 ***
## mesJune           0.206836   0.004354   47.51   <2e-16 ***
## mesJuly           0.205492   0.004270   48.12   <2e-16 ***
## mesAugust         0.235209   0.004276   55.01   <2e-16 ***
## mesSeptember      0.225107   0.004252   52.94   <2e-16 ***
## mesOctober        0.205699   0.004252   48.38   <2e-16 ***
## mesNovember       0.235935   0.004310   54.74   <2e-16 ***
## mesDecember       0.543046   0.004011  135.39   <2e-16 ***
## Feriado          -5.647972   0.059671  -94.65   <2e-16 ***
## ultimo_dia_habil  0.349682   0.003884   90.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 913454  on 1826  degrees of freedom
## Residual deviance: 101355  on 1807  degrees of freedom
## AIC: 114531
## 
## Number of Fisher Scoring iterations: 8

Veamos cómo cambia el \(R^2\):

y_pred_tr_glm<-predict(modelo_glm,type="response")
r_tr_glm<-datos_ventas_autos$und[datos_ventas_autos$Fecha<="2016-12-31"]-y_pred_tr_glm
R_tr_glm<-mean(r_tr_glm^2)
R2_tr_glm<-1-R_tr_glm/R0_tr
print(R2_tr_glm)
## [1] 0.8183819

Y en validación:

datos_val<-subset(datos_ventas_autos,subset=(Fecha>="2017-01-01" & Fecha<="2017-12-31"))
y_pred_vl_glm<-predict(modelo_glm,type="response",newdata = datos_val)
r_vl_glm<-datos_val$und-y_pred_vl_glm
R_vl_glm<-mean(r_vl_glm^2)
R2_vl_glm<-1-R_vl_glm/R0_vl
print(R2_vl_glm)
## [1] 0.6444367

FIN